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Abstract The dynamics of glacial cycles is studied in 
terms of the dynamical systems theory. We explore the 
dependence of the climate state on the phase of astro- 
nomical forcing by examining five conceptual models 
of glacial cycles proposed in the literature. The models 
can be expressed as quasipcriodically forced dynamical 
systems. It is shown that four of them exhibit a strange 
nonchaotic attractor (SNA), which is an intermediate 
regime between quasiperiodicity and chaos. Then, the 
dependence of the climate state on the phase of astro- 
nomical forcing is not given by smooth relations, but 
constitutes a geometrically strange set. Our result sug- 
gests that the dynamics of SNA is a candidate for the 
dynamics of glacial cycles, in addition to the well-known 
dynamics of quasiperiodicity and chaos. 

Keywords Glacial cycles • Astronomical theory of 
climate change • Conceptual models • Quasipcriodically 
forced dynamical systems • Strange nonchaotic 
attractors 



1 Introduction 

The climate change in the Quaternary is characterized 
by alternations between cold (glacial) and warm (in- 
terglacial) periods, which are called glacial-interglacial 
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cycles or alternatively glacial cycles. The aim of this 
article is to explore the dynamics of glacial cycles by 
examining conceptual models proposed in the litera- 
ture. 

Astronomical theories of ice ages, as typified by the 
Milankovitch theory (1941), attempt to explain glacial 
cycles based on the change in the incoming solar ra- 
diation (insolation) due to the long-term variations of 
the Earth's orbital parameters. Assuming a constant so- 
lar output and a perfectly transparent atmosphere, the 
insolation F at a given latitude and season is a func- 
tion of the astro-insolation parameters: the eccentricity 
e of the Earth's orbit, the obliquity e (the inclination of 
the equator on the ecliptic), and the climatic precession 
esina), where a) is the longitude of the perihelion mea- 
sured from the moving vernal equinox (Berger 1978). 
That is, F ~ F{e, e, esinw). 

A motion expressed by a superposition of periodic 
motions of two or more incommensurate frequencies 
is called a quasiperiodic motion. In the studies of the 
Quaternary climate changes, the astro-insolation pa- 
rameters arc often assumed to be quasiperiodic (Berger 
1978): 

e = eo + J2Ei cos(Ajt + (/)»), 

e^e* + J2A^ cos(/,i + <50, 
e sin a) = Pi sin(Q;it -I- Q). 

Then, the insolation F is also quasiperiodic in time 
(cf. Eq. (O). The insolation variation F{t) is known 
as astronomical forcing. The primary frequencies of as- 
tronomical forcing are approximately 1/19 kyr^^ and 
1/23 kyr^^ for the climatic precession, and approxi- 
mately 1/41 kyr~^ for the obliquity. The assumption 
of quasiperiodicity of astronomical forcing is consid- 
ered to be approximately valid over the last several 
million years (Berger and Loutre 1991) but may not 
be valid beyond this period because of the intrinsic 
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chaoticity of planetary motion (Lasker 1989). However, 
since our study focuses on the glacial cycles during the 
past 0.7 Myr after the mid-Pleistocene transition (Clark 
et al. 2006), we construct our theory on the assumption 
of quasiperiodic astronomical forcing. 

Hays et al. (1976) presented strong evidence for as- 
tronomical theories of ice ages. They found the pri- 
mary frequencies of astronomical forcing in the geologi- 
cal spectra of marine sediment cores. However, the dom- 
inant frequency in geological spectra is approximately 
1/100 kyr~^, although this frequency component is neg- 
ligible in the astronomical forcing. This is referred to as 
the "100 kyr problem." 

To understand the climate responses to astronom- 
ical forcing, several concepts related to the dynamical 
systems theory have been proposed in the literature. 
The linear response is a simple framework to explain the 
41, 23, and 19 kyr periodicities in the geological spec- 
tra. However, the linear response cannot appropriately 
account for the 100 kyr periodicity (Hays et al. 1976). 
Ghil (1994) explained the appearance of the 100 kyr 
periodicity as a nonlinear resonance to the combina- 
tion tone 1/109 kyr"-'^ between precessional frequencies 
1/19 kyr~^ and 1/23 kyr~^. Contrary to the linear res- 
onance, the nonlinear resonance can occur even if the 
forcing frequencies are far from the internal frequency 
of the response system. Benzi et al. (1982) proposed 
stochastic resonance as a mechanism of the 100 kyr pe- 
riodicity, where the response to a small external forcing 
is amplified by the effect of noise. See Ganopolski and 
Rahmstorf (2002) and Braun et al. (2007) for the recent 
progress related to stochastic resonance. Tziperman et 
al. (2006) proposed that the timing of deglaciations is 
set by the astronomical forcing via the phase-locking 
mechanism. If a phase locking occurs, the condition 
if = Pifi + P2f2 + ■ ■ ■ + Pn In is satisfied, where / 
is the response frequency, /i, /2, Jn are forcing fre- 
quencies, and pi, p2, Pn, and 0) arc integers 
(Romeiras et al. 1987). De Saedeleer et al. (2013) sug- 
gested generalized synchronization (GS) to describe the 
relation between the glacial cycles and the astronomi- 
cal forcing. GS means that there is a functional relation 
between the climate state and the state of astronomical 
forcing. They also showed that the functional relation 
may not be unique for a certain model. However, the 
nature of the relation remains to be elucidated. 

In this research, we analyze five conceptual models 
of glacial cycles proposed in the literature and show 
that the dependence of the climate state on the phase 
of astronomical forcing may not be given by smooth 
relations, but constitutes some geometrically strange 
set known as a strange nonchaotic attractor (SNA). 



The remainder of this article is organized as follows. 
In Section 2, we review the notions of quasiperiodically 
forced systems and attractors. In Section 3, five models 
of glacial cycles are explained. In Section 4, we analyze 
the models numerically and show possible relations be- 
tween the climate state and the phase of astronomical 
forcing. In Section 5, we mention a dynamical property 
of climate models with SNAs. Section 6 concludes this 
article. 

2 Quasiperiodically forced systems and 
attractors 

2.1 Quasiperiodically forced systems 

We consider the dynamics of glacial cycles in the frame- 
work of quasiperiodically forced dynamical systems. Let 
= M^/(27rZ)^ be an TV-dimensional torus. In gen- 
eral, quasiperiodically forced dynamical systems can be 
expressed as 

(1) 

jc = g(x,6i), xeM^^ 

where 6 — {61,62, ■■■,0n)'^ is the phase of the drive 
subsystem, x = {xi,X2, ■■■,xm)'^ is the state of the re- 
sponse subsystem. g(x, 6) a periodic function in each 
phase 6i, and uj = (wi, 0^2, '^tv)'^ is a vector of in- 
commensurate frequencies such that kiui + k2U)2 4- ■ • ■ -I- 
k^oJiq = does not hold for any set of integers, fci, fc2, 
. . . , kiq , except for the trivial solution ki ~ k2 — ■ ■ ■ = 
k^ = 0. In the astronomical theory of climate change, 6 
corresponds to the phase of astronomical forcing F{t), 
and x corresponds to the climate state, as described in 
Section 3. 

2.2 Attractors 

Like most conceptual models of glacial cycles, we as- 
sume the existence of attractors in systems ([T|). An 
attractor is a compact set with a neighborhood such 
that, for almost every initial condition in this neighbor- 
hood, the limit set of the orbit as time tends to infinity 
is the attractor (Grebogi et al. 1984). The closure of 
set of initial conditions which approach an attractor as 
time tends to inifinity is called the basin of attraction 
of the attractor. Note that the attractor of system ^ 
is a subset of x M*^. Thus, the dependence of the 
climate state on the phase of astronomical forcing is 
represented by this attractor, after a certain transient 
time. In quasiperiodically forced systems ((!]) , the follow- 
ing types of attractors typically appear: N -frequency 



Dynamics between order and chaos in conceptual models of glacial cycles 



3 



quasiperiodic attractors, L-frequency quasiperiodic at- 
tractors with L > N, Strange nonchaotic attractors 
(SNAs) and chaotic attractors. 

Lyapunov exponents give the mean exponential rates 
of divergence of nearby orbits. Although the Lyapunov 
exponents depend on initial conditions, they have the 
same set of values for almost every initial condition, 
under general circumstances. Such initial conditions or 
orbits which give the same set of values of the Lya- 
punov exponents are called typical. A chaotic attractor 
is an attractor for which typical orbits on the attrac- 
tor have a positive Lyapunov exponent (Grebogi et al. 
1984). That is, almost all orbits on the attractor have 
sensitive dependence on initial conditions. In quasiperi- 
odically forced systems ([T]), an attractor is chaotic if 
the maximal conditional Lyapunov exponent (CLE) A 
is positive for typical orbits on the attractor, and non- 
chaotic otherwise. The maximal CLE A is defined as 

A = lim - In ] ^ f „| \ , 
t^oct |(5x(0)r 

where (5x(t) is an infinitesimal displacement of the or- 
bit from x(i). The evolution of the infinitesimal dis- 
placement is given by the variational equation (5x — 
DxG(x, e)6K. 

An TV-frequency quasiperiodic attractor is an at- 
tracting iV-dimensional torus on which the motion is 
quasiperiodic. A^-frequency quasiperiodic attractors are 
characterized by a negative maximal CLE, A < 0. 

An L-frequency quasiperiodic attractor with L > N 
is an attracting L-dimensional torus on which the mo- 
tion is quasiperiodic. L-frequency quasiperiodic attrac- 
tors with L > N are characterized by a maximal CLE 
of zero, A = 0. 

A strange nonchaotic attractor (SNA) is a strange 
attractor for which typical orbits have nonpositive Lya- 
punov exponents (Grebogi et al. 1984). The strange at- 
tractor is defined as an attractor which is not a finite set 
of points and is neither piecewise differentiable curve or 
surface, nor a volume bounded by a piecewise differen- 
tiable closed surface. SNAs are characterized by a nega- 
tive maximal CLE, A < 0. It is proven that some SNAs 
are represented by a closure of graphs of almost every- 
where discontinuous functions (Keller 1996; Alseda and 
Costa 2009). 

SNAs lie in an intermediate regime between quasiperi- 
odicity and chaos. The power spectrum of an SNA is 
intermediate one between quasiperiodicity and chaos, 
namely, a singular continuous (fractal) spectrum (Feudel 
et al. 1996). SNAs are multifractal objects in the sense 
that the capacity dimension is larger than the informa- 
tion dimension (Ding et al. 1989; Hunt and Ott 2001). 
For comprehensive reviews of SNAs, refer to Prasad et 
al. (2001) as weh as Feudel et al. (2006). 



In the remainder of this article, we will not en- 
counter examples of L-frequency quasiperiodic attrac- 
tors with L > N, and sometimes A^-frequency quasiperi- 
odic attractors are simply referred to as quasiperiodic 
attractors. 



2.3 Examples of A^-frequency quasiperiodic attractors 
and an SNA 

A quasiperiodically forced map is regarded as a stro- 
boscopic map obtained by sampling the state of the 
continuous-time system ([T]) at time intervals correspond- 
ing to one of the forcing frequencies. Let us consider a 
quasiperiodically forced map studied by Glendinning 
(2002): 



0{n + I) = e{n) + LO (modi), 

x{n + 1) = 2a{a + cos27r6'(n)) tanh.T(n), 



(2) 



where {9, x) e [0,1) x R, n e Z, uj e R \ Q, and a and 
a arc non-negative parameters. 

Figure 1(a) shows the quasiperiodic attractors for 
cr = 1.2, a = 1.001, and uj = (V^ - l)/2 (solid lines). 
Each attractor is represented by a graph of a smooth 
function. If we denote the attractor in the region a; > 
hy X = (f>{0), the attractor in the other region is x = 
^<j){0) because of the symmetry x — )■ —x. The line a; = 
is an unstable invariant set, namely, a repellor. In this 
case, the repellor constitutes the basin boundary of two 
attractors. 

Figure 1(b) shows the SNA for a = 1.2, a = 0.8, 
and uj = — l)/2. The SNA is represented by the 
closure of two graphs x = (p{9) and x = —(p{9) having 
following properties (Glendinning 2002; Keller 1996): 

(i) if {9) > for ah 9 e [0,1); 

(ii) the union of two graphs x = f{9) and x = —ip{9) 
is invariant under the map; 

(iii) typical conditional Lyapunov exponents on the graphs 
are negative; 

(iv) (p{9) is discontinuous at almost all 9 E [0, 1); 

(v) (p{9) is upper semi-continuous; 

(vi) ip{9) = on a dense, but measure zero, set of values 
of 6* G [0,1). 

Unlike quasiperiodic attractors, SNAs contain some 
repellors in themselves (Stark 1997). See the repellor 
a; = embedded in the SNA in Fig. [ijb). Due to such 
repellors, the orbits on SNAs become temporarily un- 
stable and sensitive to perturbations, as shown in Sec- 
tion 5. 
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Fig. 1 Attractors in map a Quasiperiodic attractors for 
a = 1.2, a = 1.001, and uj = {\/5—l)/2 (upper and lower solid 
lines)DThe maximal CLE is A = —0.3677 for each attractor. 
b SNA for cr = 1.2, a = 0.8, and lj = {V5 - l)/2 (dots). The 
maximal CLE is A = —0.2124. In each panel, the dashed line 
a; = is the unstable invariant set 



2.4 Numerical characterization of strangeness of 
attractors 

SNAs and A^-frequency quasiperiodic attractors can be 
numerically distinguished by the phase sensitivity ex- 
ponent (Pikovsky and Feudel 1995) or by the parame- 
ter sensitivity exponent (Nishikawa and Kaneko 1996). 
First, let us consider the case of the continuous-time 
system ([T]). The phase sensitivity function r{t) of a 
variable Xj with respect to a variable 9i is defined as 



mm 

c(o),e(o)} 



max 

0<T<t 



The derivative dxj{T)/d9i{0) characterizes the sensitiv- 
ity of Xj with respect to a change of 9i . This derivative 
is obtained by integrating 



d dxj 
di 89,(0) 



E 

k=l 



dgj{x,9) dxk dg, (x, ( 



dxk d9,{0) d9,{0) 



(3) 



along the solution {9{t)^ x(t)). For SNAs, the phase sen- 
sitivity function grows as r{t) ~ with > 0, where 
^ is called the phase sensitivity exponent. On the other 
hand, for A^-frequency quasiperiodic attractors, r{t) is 
bounded and ^ = 0. 

Assume that function g(x, 9) contains some param- 
eter, say a. The parameter sensitivity function ^{t) of 
variable Xi with respect to parameter a is defined as 

dXj{T) 



9(t) = min < max 

{x(0),6l(0)} I 0<r<t 



da 



The derivative dxj (r) / da characterizes the sensitivity 
of Xj with respect to a change of a. This derivative is 
obtained by integrating 



d dx j 



M 

dt da ^ 



dgj{x,9) dxk , 9gj(x, 



dxk 



da 



da 



(4) 



along the solution {9{t), x(t)). For SNAs, the phase sen- 
sitivity function grows as fl'(t) ~ with v > 0, where 
V is called the parameter sensitivity exponent. On the 
other hand, for A^-frequency quasiperiodic attractors, 
^{t) is bounded and u ~0. 

For discrete-time systems, such as Eq. the meth- 
ods for characterizing the strangeness of attractors are 
almost the same except that time evolutions correspond- 
ing to Eqs. ^ and Q become time-discrete (Pikovsky 
and Fcudel 1995; Nishikawa and Kaneko 1996). 

Figurc[2]shows the behaviors of the phase sensitivity 
function r[n) and the parameter sensitivity function 
^{n) with respect to a for the attractors in Fig. [1] The 
functions r{n) and S'(n) saturate for large n for the 
quasiperiodic attractors, and on the other hand, show 
power- law divergences with exponents pi = 1.03 and 
— 0.99, respectively, for the SNA. 

3 Models 

3.1 Astronomical insolation forcing 

Milankovitch (1941) considered the high-latitude north- 
ern hemisphere summer insolation as a decisive factor 
for glaciations. We use the insolation F{t) at 65°N on 
the day of the summer solstice in the simulations of 
all the models. For simplicity in computation, the in- 
solation F{t) is calculated using the formula given in 
De Saedeleer et al. (20130: 

1 ^ 

F{t) = - ^[si sm{LL!,t) + Ci cos{uj,t)], 

i=l 

^ According to De Saedeleer et al. (2013), the formula is 
valid over the past 1 Myr, and its mean error compared to 
Lasker et al. (2004) is 6.7 W/m^ with peaks at 27.5 W/m^. 
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Table 1 Sets of parameter values used for the SM90 and 
MS90 models 



1 2 3 4 5 6 7 



a=1.001 
a=0.8 



(b) 



-2 



3 4 5 



Fig. 2 a Phase sensitivity function r(n) for the quasiperi- 
odic attractors of map with a = 1.2, a = 1.001, and 
LJ = (\/5— l)/2 {open square) and for the SNA with a = 1.2, 
Of = 0.8, and cj = (\/5 — 1)/2 [filled circle), b Phase sensitivity 
function 'I' in) for the same attractors 



where iV = 35 is the number of different frequencies LOi 
included in the sum. The values of frequencies uJi are 
those in Berger (1978). The coefficients Si and ci are 
extracted from the solution of Berger (1978) by per- 
forming a linear regression (De Saedeleer et al. 2013). 
To make the paper self-contained, the values of Ui, Si, 
and Ci are listed in Table [3] in Appendix 1. The pa- 
rameter a is a scaling constant, the value of which is 
different among models. The top panel of Fig. 0] shows 
the insolation F{t) normalized to have a zero-mean and 
unit variance (with a = 23.58). 

Introducing variables 9i{t) = ujit [i ~ 1,2,...,A^), 
the above equation for F{t) can be written as 



Set 


P 


q 


r 


s 


u 




w 


A 


1.0 


2.5 


0.9 


1.0 


0.6 


0.2 


0.5 


B 


0.82 


4.24 


0.95 


0.53 


0.32 


0.02 


0.66 


C 


1.0 


1.2 


0.8 


0.8 


0.7 


0.0 


0.0 



3.2 Conceptual models of glacial cycles 
3.2.1 Saltzman-Maasch model 

Saltzman and Maasch (1990) introduced a global model 
of the late Cenozoic climate. We will refer to this model 
as SM90. The model consists of differential equations 
for the global ice volume x, the atmospheric CO2 con- 
centration y, and the global mean deep water tempera- 
ture z, where x, y, and z are dimensionless and denote 
deviations from 1 Myr averages. Neglecting stochastic 
perturbations, the model is written in the following de- 
terministic from: 



Tx ~ —X — y — vz ~ uF{t), 
.2 



ry ~ —pz + ry + sz 
Tz = —q{x + z) 



wyz 



(6) 



where t is the time in kiloyears, and r = 10 kyr. The 
parameters p, g, r, s, and w are freely chosen, but 
the remaining parameters u and v are functions of the 
other parameters. The insolation forcing Fit) is given 
by Eq. © with a = 23.58. 

Saltzman and Maasch (1990) simulated glacial cy- 
cles for the last 500 kyr by choosing the parameter val- 
ues as set A in Table 1. Hargreaves and Annan (2002) 
estimated the parameter values of SM90 as set B using 
a data assimilation technique. 

3.2.2 Maasch-Saltzman model 

Prior to the SM90 model, Maasch and Sahzman (1990) 
studied one of their previous models. This older model, 
which we call the MS90 model, can be obtained from the 
SM90 model by setting v = w ^ 0. The other parameter 
values are chosen as set C in Table [TJ 



N 



F{t) = i V[s, sin 6*,; (t) + c, cos 61,(0], 
a ^-^ 

1=1 

&»(0) = (i = l,2,...,iV). 



(5) 



Note that Qi arc the phase variables of astronomical 
forcing in climate models. 



3.2.3 Crucifix-De Saedeleer-Wieczorek model 

One of the simplest systems that exhibit self-sustained 
oscillations is the van der Pol oscillator. With slight 
modifications and addition of an insolation forcing to 
the van der Pol oscillator. Crucifix (2012) introduced 
the following model: 



Tx^-[y + P + jF{t)], 
ry = -a[y^/3 - y - x], 



(7) 
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where x denotes the global ice volume, and y is a con- 
ceptual variable for realizing self-sustained oscillations. 
The insolation forcing F{t) is given by Eq. ([5|) with 
a = 11.77. For the unforced case 7 = 0, Eq. ([7]) has one 
stable equilibrium state for > 1 and a limit cycle 
oscillation for |/?| < 1. Following Crucifix (2012), the 
parameters are set as a = 30.0, /3 ~ 0.75, 7 = 0.4, and 
T = 36.0 kyr. 

The dynamical properties of Eq. ([7]) were exten- 
sively studied by De Saeldeleer, Crucifix, and Wieczorek 
(2013). Thus, we refer to Eq. © as the CSW model. 

3.2.4 Paillard-Parrenin model 

Paillard and Parrenin (2004) proposed an oscillator model 
of glacial cycles involving the Antarctic ice-sheet influ- 
ence on bottom water formation. The model, which we 
call the PP04 model, consists of equations for the global 
ice volume x, extent of the Antarctic ice sheet y, and 
atmospheric CO2 concentration z: 

TxX = —pz — qF{t) + r — X, 
Tyy = x- y, 



TzZ = aF{t) - Px + '-fH{-w) + 5- z, 
w = ax ~ by — cl(t) + d, 



(8) 



where F{t) is the summer solstice insolation at 65°N 
given by Eq. ([S|) with a = 23.58, and I{t) is the in- 
solation at 60° S on February 21. The Heaviside func- 
tion in the original PP04 model is approximated as 
H{x) — 1/2 -f arctan(1000x)/7r for simplicity of the 
dynamical systems analysis. The parameters are set as 
= 15 kyr, Ty = 12 kyr, = 5 kyr, p = 1.3, q = 0.5, 
r = 0.8, a = 0.15, 13 = 0.5, 7 = 0.5, 6 = 0.4, a = 0.3, 
b = 0.7, c = 0, and d ~ 0.27. The choice c = does not 
change the timing of terminations; it only affects the 
qualitative agreement with geological records (Paillard 
and Parrenin 2004) . 

3.2.5 Imbrie-Imbrie model 

Imbrie and Imbrie (1980) proposed a piecewise-linear 
model for the glacial cycles: 

f {l + b)(F{t)-x) if F(t) >x, 
\{l^b){Fit)-x) ffF(t) <x, 

where x corresponds to a loss in the global ice volume. 
The variable F{t) is assumed to be the summer solstice 
insolation at 65°N given by Eq. ([S|) with a — 23.58 (cf. 
Paillard 2001). The parameters are tuned as r = 17 kyr 
and b = 0.6. For simplicity of the dynamical systems 
analysis, we deal with a smoothed system. 



If the steepness parameter k is large enough, the solu- 
tion of Eq. ^ converges to that of the original model. 
We set fc = 10 in this study. 



4 Dynamical systems analysis 

4.1 Synchronization under the same astronomical 
forcing 

In what follows, x{t) represents the climate state at 
time t for each model described in the previous sectiorl^. 
First, we check the dependence of solutions on initial cli- 
mate states under the same astronomical forcing F(t). 
We simulate lO'' solutions x(*)(i) {i = 1, 2, 10-*), 
starting from random initial states x^*^(to) at time to = 
—20 Myr, under the same astronomical forcing F{t). 
For these solutions, we calculate the maximum distance 
from the average of the ensemble. 



At) ~ max 

l<i<10* 



(0-(xW(i)) 



TX 



[1 + 6tanhfc(F(t) - x)]{F{t) - x). 



(9) 



We assume that typical solutions converge to the same 
solution, i.e., synchronize, under the same astronomical 
forcing if A^ax{t) — >■ as t — )• 00. Figure |3] shows the be- 
havior of Z\„iax(i) for each model. In each model except 
for the SM90-B model, typical solutions synchronize af- 
ter a certain transient time. This result implies that 
the attractor is unique at least except for the SM90-B 
model. It should be mentioned that the transient time 
is longer than the duration of the Quaternary, approxi- 
mately 2.58 Myr, in the SM90-A and MS90 models (see 
Fig.lll). 



4.2 Simulated ice volume over the last 700 kyr 

Figure HI shows the simulated normalized ice volume 
V{t) over the last 700 kyr (solid line) for each model. 
The S^^O record of the LR04 stack is also shown by a 
dashed line for comparison (Lisiecki and Raymo 2005). 
In the SM90-A, PP04, and Imbrie models, synchronized 
solutions after each transient period are well accorded 
with the 6^^0 record; however, in the MS90 and CSW 
models, synchronized solutions after each transient pe- 
riod are not accorded with the S^^O record, but only 
some transient solutions agree with the 6^^0 record. 
In the SM90-B model, solutions starting from different 
initial states do not synchronize. This result indicates 
that the SM90-B model exhibits chaotic solutions. 

Actually we find coexistence of a quasiperiodic at- 
tractor and a chaotic attractor for the SM90-B model. 

^ For example, x = {x, y, z)"^ for the SM90 model, and 
X = {x)'^ for the Imbrie model. 
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Fig. 3 Convergence of 10* solutions for each model: 
The maximum distance from the average of the ensemble 
^max(t) = maxi<i<io4 |x(*^(t) — (x'*'(t))[ vs. elapsed time 
t — Iq. All solutions start from random initial states at time 
to = — 20 Myr. Note that the time scale is different for each 
panel 



Fig. 4 Top panel: Normalized insolation F{t) at 65° N on the 
day of summer solstice. Other panels: Simulated normalized 
ice volume V{t) {solid lines) for the models described in Sec- 
tion 3. The S^^O record of the LR04 stack is also shown for 
comparison {dashed lines). Post-transient solutions are con- 
sidered for the SM90-A, SM90-B, PP04, and Imbrie models, 
and transient solutions are considered for the MS90 and CSW 
models 



The quasiperiodic attractor originates in a stable equi- 
librium of the unforced system (x, y, z) =(-0.282, 0.277, 
0.282). The basin of attraction of this quasiperiodic at- 
tractor is much smaller than that of the chaotic attrac- 
tor, and the amplitude of the quasiperiodic solution is 
also much smaller than that of the chaotic solutions. 
Thus, we do not focus on this quasiperiodic solution. 



4.3 Types of attractors 



Table 2 Characteristic exponents and types of attractors 
for the different models. The asterisk (*) indicates that the 
sensitivity exponents are undefined for the SM90-B model. 



Models 


A [kyr-i] 


At 


V 


Attractor 


SM90-A 


-0.0014 


1.06 


1.15 


SNA 


SM90-B 


-h0.0020 


* 


* 


Chaotic 


MS90 


-0.0005 


0.95 


1.09 


SNA 


CSW 


-0.016 


1.10 


1.13 


SNA 


PP04 


-0.0093 


1.09 


1.08 


SNA 


Imbrie 


-0.046 


0.0 


0.0 


Quasiperiodic 



We analyze types of attractors using the maximal CLE 
A, the phase sensitivity exponent ^, and the parame- 
ter sensitivity exponent v. As shown in Tabic 2, the 
maximal CLE A is positive for the SM90-B model and 
is negative for all the other models. Thus, the SM90- 
B model is chaotic, and all the other models are non- 



chaotic. This result is consistent with that of conver- 
gence experiment. 

The next question is whether the attractors are strange 
or not. Figurc[51^a) shows the phase sensitivity functions 
r(i) of X with respect to a phase variable for the 
models with nonchaotic attractors, and Fig.[SJb) shows 
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Fig. 5 a Phase sensitivity function r{t) of x with respect 
to a phase variable 84. b Parameter sensitivity function 'Z'(t) 
of X with respect to u for the SM90-A and MS90 models, 
parameter /3 for the CSW model, and parameter b for the 
Imbrie model 



the parameter sensitivity functions S'(t) of x with re- 
spect to parameter u for the SM90-A and MS90 models, 
parameter /3 for the CSW model, and parameter b for 
the Imbrie model. Here, the functions r{t) and tf'(t) are 
calculated for a set of 20 initial conditions. For the Im- 
brie model, the functions r{t) and 'F{t) saturate as time 
increases, but for the other models, the functions r{t) 
and ^(t) grow in power-laws with exponents /i > and 
v > 0, respectively. Table 2 shows the values of 11 and 
u, which are estimated from the power-law behaviors 
for t > 10''. Therefore, the attractors of the SM90-A, 
MS90, CSW, and PP04 models are SNAs, and the at- 
tractor of the Imbrie model is a quasiperiodic attractor. 



attractors in climate models have at least 35 dimen- 
sions, which are too high-dimensional to be visualized. 
To capture the geometric features of these attractors, 
we illustrate approximated attractors by a truncation 
of the astronomical forcing. Leaving only the 23 kyr 
(i = 1) and 19 kyr {i = 3) precession terms and the 
41 kyr (i = 4) obliquity term, we simplify the astro- 
nomical forcing F{t) as 



1 



0(1,3,4} 



[si sin(a;ii) 4- Ci cos(a'it)] . 



i6{l,3,4} 



These three frequency components constitute 78% of 
the original astronomical forcing with respect to stan- 
dard deviations. The scaling constant is a{i^3^4} — 9.08 
for the CSW model and a^i 3 4} = 18.3 for the other 
models. 

In spite of the simplification of the astronomical 
forcing, we obtain the same type of attractor for each 
model although the values of characteristic exponents 
are slightly changed: an SNA with A ~ —0.0004, n = 
1.39, and z/ = 1.43 for the SM90-A model, a chaotic 
attractor with A « 0.0025 for the SM90-B model, an 
SNA with A « -0.00012, ^ = 1.09, and i' = 1.35 
for the MS90 model, an SNA with A « -0.0088, fi = 
1.13, and i' = 1.10 for the CSW model, an SNA with 
A « -0.0089, n = 1.07, and z/ = 0.94 for the PP04 
model, and a quasiperiodic attractor with A « —0.046 
and ve ^ = for the Imbrie model. 

The stroboscopic map is useful for visualizing at- 
tractors having three or more degrees of freedom. Ob- 
serving the system at the moments t„ = 27rn/a;i {n G 
Z), where a phase Oi{t) attains a constant value, the 
system reduces to a stroboscopic map {03{tn), Oi{tn), 
^{tn)) ^-> {O^itn+i), 94{tn+i), x(t„+i)). Figure H] shows 
the attractors of the simplified models. Each panel presents 
the projection of the attractor on the space (6'3(t„)/27r, 
e4{tn)/2TT, x{tn)). The SM90-A, MS90, CSW, and PP04 
models exhibit SNAs, and the dependence of the cli- 
mate state on the phase of astronomical forcing is not 
given by smooth relations, but constitutes a geometri- 
cally strange set. Moreover, the dependence seems to be 
discontinuous almost everywhere. For the Imbrie model, 
the attractor is quasiperiodic, and the dependence is 
given by a smooth graph. For the SM90-B model, the 
attractor is chaotic, and the dependence constitutes a 
geometrically strange set, as in the case of SNAs. 



4.4 Illustration of attractors by simplification of 
astronomical forcing 

Since the astronomical forcing F{t) in Eq. ([5]) has high 
degrees of freedom, namely, = 35 in our study, the 



5 A dynamical property of SNAs: Sensitivity to 
perturbations 

If a climate model written as Eq. ([1]) has a chaotic at- 
tractor, solutions sensitively depend on initial climate 
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Fig. 6 Attractors of each model for a simplified astronomical forcing. Each panel shows the projection of the attractor on the 
space 04/27r, x) 



states x(io). On the other hand, if a chmate model has 
an SNA or a quasiperiodic attractor, we obtain one or 
several synchronized solutions after a certain transient 
time. Then what is the difference between SNAs and 
quasiperiodic attractors? 

Khovanov et al. (2000) stated that SNAs and quasiperi- 
odic attractors differ in sensitivity to perturbations. We 
illustrate the difference by simulating the CSW and Im- 
brie models in the presence of noise. In the former, we 
add independent white Gaussian noises ^i{t) [i = x 
or y) to the right-hand sides of Eq. d?]), respectively, 
where {i,{t)) = 0, {i,{t)ij{t')) = a^5,,5{t - t'), and 
a = 0.1 kyr ' . In the latter, we add a white Gaussian 
noise ^xit) with the same average and standard devia- 
tion to the right-hand side of Eq. 

Figure[71^a) shows two trajectories of the CSW model 
simulated with and without noises, and Fig.[7Jd) shows 
those of the Imbrie model. The CSW model is more 
sensitive to noise than the Imbrie model because of the 
local instability of SNA. To characterize the local in- 
stability of attractors, let us define a finite-time CLE 



1 \6^{t + T)\ 



The finite-time CLE Xxit) gives the rate of exponen- 
tial divergence of nearby orbits during the time inter- 
val [t, t + T]. Note that Xrit) converges to the maximal 
CLE A as T -> oo. For SNAs, the finite-time CLE Xrit) 



can be positive even for large t with nonzero proba- 
bility (Pikovsky and Feudel 1995). Thus, the orbits of 
SNAs are sensitive to perturbations in time intervals 
with At(0 > 0. Figures [71[a)-(c) show that the per- 
turbed trajectory occasionally deviates after the finite- 
time CLE Xxit) becomes positive, where T ~ 70 kyr. 
This local instability may be attributed to a repellor 
embedded in SNA as mentioned in Section 2. Such a 
temporary loss of synchronization has been shown also 
in another model of glacial cycles (Crucifix 2011). 

On the other hand, for quasiperiodic attractors, the 
finite-time CLE Xrit) takes only negative values for 
large T. As a result, the orbits of quasiperiodic attrac- 
tors arc robust to perturbations in comparison with 
those of SNAs. Figures [71[d)-(f) show that the Imbrie 
model is robust to the same magnitude of noise that 
can cause deviations of trajectories in the CSW model. 



6 Conclusion 

We analyzed five conceptual models of glacial cycles 
in the framework of the dynamical systems theory. The 
SM90-A, MS90, CSW, and PP04 models exhibit a strange 
nonchaotic attractor (SNA), and the dependence of the 
climate state on the phase of astronomical forcing is 
not given by smooth relations, but constitutes a geo- 
metrically strange set. The SM90-B model exhibits a 
chaotic attractor, and the dependence of the climate 
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Fig. 7 Sensitivity of climate models to noise: a Ice volume 
x(t) simulated by the CSW model with noise (black) and 
without noise {red), b Distance A[t) in the space {x,y) be- 
tween perturbed and unperturbed trajectories corresponding 
to panel a. c The finite-time CLE Xxit) of the unperturbed 
trajectory in panel a. T = 70 kyr. d Loss of ice volume x{t) 
simulated by the Imbrie model with noise {black) and without 
noise {red), e Distance A{t) between perturbed and unper- 
turbed trajectories of panel d. f The finite-time CLE XT{t) 
of the unperturbed trajectory in panel d. T = 70 kyr 



forcing. In their article, the inhibition of chaotic behav- 
ior refers to the appearance of statistical periodicity for 
the ensemble of orbits, and it is different from the non- 
chaoticity based on nonpositive Lyapunov exponents in 
this study. 

The analysis in this paper was restricted to spatially 
zero-dimensional models that can be approximated by 
smooth dynamical systems. A hybrid dynamical sys- 
tem is a system that has both continuous and discrete 
dynamics (Aihara and Suzuki 2010). A number of ice 
age models are expressed as hybrid dynamical systems, 
more precisely, hybrid automat^ (van der Schaft and 
Schumacher 2000). For example, the models by Paillard 
(1998), Parrenin and Paillard (2003), and Ashkenazy 
and Tziperman (2004) are hybrid automata. The Im- 
brie and PP04 models are also hybrid dynamical sys- 
tems but can be approximated by smooth dynamical 
systems because they are especially piecewise ajfine sys- 
tems. The description of the dynamics appearing in hy- 
brid automata models or in more complicated models 
with spacial dimensions is still unclear. This problem 
will be studied in the future. 
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Appendix 1: Parameter values for the insolation 
formula 



state on the phase of astronomical forcing is not given 
by smooth relations, but constitutes a geometrically 
strange set. The Imbrie model exhibits a quasiperiodic 
attractor, and the dependence of the climate state on 
the phase of astronomical forcing is given by a smooth 
graph. These results suggest that the dynamics of SNAs 
is a candidate for the dynamics of glacial cycles, in ad- 
dition to the well-known dynamics of quasiperiodicity 
and chao^. Since SNAs have repellors in themselves, 
the orbits on SNAs are locally unstable. Therefore, if 
the real climate system has an SNA, the orbit will be 
temporarily sensitive to perturbations when the finite- 
time conditional Lyapunov exponent is positive. 

Brindlcy and Kapitaniak (1992) specified an inhibi- 
tion of chaotic behavior due to the quasiperiodicity of 

^ Chaos in conceptual models of glacial cycles has often 
been mentioned in the literature (e.g. Nicolis 1987; Saltzman 
and Verbitsky 1992, 1993; Ghil 1994; Huybers 2009). 



The parameter values for the insolation formula (O 
given by De Saedeleer et al. (2013) are listed in Table[3] 
to make the paper self-contained. The frequencies and 
coefficients are arranged in descending order of their 



power sf -\- cf. 
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